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,-0 Summary. We describe a new method to include magnetic fields into smooth par- 

ticle hydrodynamics. The derivation of the self-gravitating hydrodynamics equations 
from a variational principle is discussed in some detail. The non-dissipative magnetic 
field evolution is instantiated by advecting so-called Euler potentials. This approach 
enforces the crucial V B = O-constraint by construction. These recent developments 
are implemented in our three-dimensional, self-gravitating magnetohydrodynamics 
code MAGMA. A suite of tests is presented that demonstrates the superiority of 
this new approach in comparison to previous implementations. 
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1 Specific astrophysical requirements 

Astrophysical simulations have their specific requirements which differ in 
CD many respects from those of other branches of computationally intense re- 

search fields. The dynamics of self-gravitating gas masses plays a prominent 
role throughout astrophysics, but it is usually only one of several ingredients 
and it is often necessary to account for additional physical processes such 
as radiative transfer, nuclear burning or the evolution of magnetic fields to 
J> address questions of astrophysical interest. These additional processes often 

involve intrinsic length and time scales that are dramatically different from 
those of the gas dynamical processes making astrophysical problems prime 
examples of multi-scale and multi-physics challenges. 

Since fixed boundaries are usually absent, flow geometries are determined by 
the interplay between different physical processes such as gas dynamics and 
(self-) gravity which often leads to complicated, dynamically changing flow ge- 
ometries. Therefore, many problems require flexible numerical schemes such 
as adaptive mesh refinement or completely meshfree, Lagrangian methods. 
Each of these methods has its stengths and weaknesses and the choice of the 
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Figure 1. Two snapshots from the simulation of the tidal disruption of a white 
dwarf by a black hole. 



best-suited method can usually save a tremendous amount of effort. 
An astrophysical example of such an intrinsic multi-scale and multi-physics 
problem is shown in Fig. I 3 . It shows two snapshots from a simulation of the 
tidal disruption of a white dwarf star by a black hole. The initially spherical 
star becomes strongly distorted while passing the black hole (left panel), it 
is heavily compressed compressed and shock-heated which triggers very rapid 
nuclear reactions whose energy release leads to the thermonuclear explosion of 

3 Astrophysical implications of this topic are discussed in [1, 2], details of the nu- 
merics can be found in [3] 
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the white dwarf. In order to follow this process for each of the computational 
fluid particles a nuclear network [4] is evolved on-the-fly together with the 
hydrodynamics . 

In many astrophysical problems the numerical conservation of physically con- 
served quantities determines the success and the reliability of a numerical 
simulation. Consider, for example, a molecular gas cloud that collapses under 
the influence of its own gravity to form stars. If in the simulation angular mo- 
mentum is artificially dissipated, say due to too coarse a mesh discretization, a 
collapsing, self-gravitating portion of gas may form just a single stellar object 
instead of a multiple system of stars and it will thus produce a qualitatively 
wrong result. The "exact" 4 conservation of mass, energy, linear and angu- 
lar momentum is -besides its natural adaptivity- one of the main strengths 
of the smoothed particle hydrodynamics (SPH) method. This exact conser- 
vation can be "hardwired" into SPH's evolution equations via symmetries in 
the fluid particle indices. Originally this was done -successfully, but somewhat 
arbitrarily- by hand [5, 6, 7], but more recently it was shown [8, 9, 10, 11] how 
the correct symmetries follow elegantly and stringently from a discretized fluid 
Lagrangian and the Euler-Lagrange variational principle. 
In the following, we will review the derivation of the self-gravitating SPH 
equations from a Lagrangian, see Sect. 2.1. We will also discuss in detail how 
to implement magnetic fields via so-called Euler potentials, see Sect. 2.2. This 
approach is similar to evolving a vector potential and enforces the crucial 
V • B = 0-constraint which otherwise poses a severe challenge for particle 
methods. These new dcveloments arc implemented in our self- gravitating, 
three-dimensional magnetohydrodynamics code MAGMA, which is described 
in Sect. 3. 



2 Guiding principles 

2.1 Ideal smoothed particle hydrodynamics (SPH) 

The smoothed particle hydrodynamics method (SPH) had originally been de- 
veloped in the astrophysical context to simulate the formation of stellar binary 
systems via fission[12] and the structure of non-sperical stars[13]. While the 
initial 3D simulations used 80 (Gingold and Monaghan) and 100 SPH parti- 
cles (Lucy) today's state of the art cosmological SPH simulations have reached 
particle numbers in excess 10 9 , see e.g. [14]. This is only in part due to the in- 
crease in hardware performance, also the simulation techniques (in particular 
the treatment of self-gravity) have become continuously more sophisticated 
and much effort has been invested to parallelize 3D codes on various comput- 
ing platforms. Also the formulation of the SPH equations has come a very long 

4 "Exact" means up to possible effects from the numerical integration of the re- 
sulting ODEs or from using approximative forces, say from a tree or some other 
Poisson-solver. 
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way from the initial straight-forward discretisation of the Lagrangian gas dy- 
namics equations to its most recent formulation that follows stringently from 
a discretized ideal fluid Lagrangian. 

Here we will give a brief overview over an older SPH- formulation, but we 
will mainly focus on an approach that is based on a derivation from a discre- 
tised Lagrangian. This latter approach naturally introduces so-called "grad-h" 
terms that result from changes in the smoothing lengths of the SPH particles. 

"Vanilla Ice" SPH 

The approximation of function values and derivatives via a kernel summation 
is at the heart of SPH. If the values of a function / are known at a set 
of discrete points ( "particles" ) labelled by b, the SPH approximation of the 
function / at position r is given by [6, 7, 15] 5 

f{r) = Y,—hW{r-r b ,h), (1) 
b pb 

where mt is the (usually constant) particle mass, p b is the mass density and 
W is a kernel function whose width is determined by the smoothing length h. 
Essentially all astrophysical SPH codes use the cubic spline kernel suggested 
in [16]. Kernel functions with compact support are preferable since they re- 
strict the SPH-summations to a local set of neighbours. For the conservation 
properties it is convenient to have "radial" kernels, 

W(r-r b ,h) = W(\\r-r b \\,h), (2) 

so that 

VaW bk = V b W kb (8 ba - 4a), (3) 

and 



V a W ab = — e ab , (4 

dr ab 

\\r ab \\, W ab = W(r ab ,h) and e ab = r ab /r ab . This 
V a W ab = -X7 b W ab (5) 



where r ab = r a - r b , r ab = 
immediately leads to 

and 



- = v ab ■ V a W ab , 



(6) 



5 Note that we do not specify at this point which h is used. For this "vanilla ice" 
SPH the h that enters the kernel should be a symmetric combination of the 
smoothing lengths of the involved particles. This will be explained in more detail 
below. For simplicity, we are omitting the subscript h in what follows. We also 
drop the distinction between the function to be interpolated and the interpolant, 
i.e. we use the same symbol / on both sides of the following equation. 
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with v a b = v a — Vb being the velocity difference between particle a and b. 
Eq. (1) can be applied in particular to the mass density itself which then reads 

p(r) = ^2m b W(r-r b ,h). (7) 

b 

The gradient of a function is approximated in SPH by taking the exact deriva- 
tive of the approximant: 

V/(r) = V— AV^r-r^). (8) 
b P b 

The most straightforward and historically first taken approch is to apply this 
set of rules to the Lagrangian form of the ideal hydrodynamics equations: 

§--*•., (9) 

which express the conservation of mass, momentum and energy. Here, P is 
the thermodynamic pressure, / abbreviates body forces and u is the thermal 
energy per mass. 

To briefly illustrate the dependence of conservation on the symmetry of the 
particle indices let us apply Eq. (8) straightforward to the pressure gradient 
in Eq. (10) (and assume vanishing body forces) to obtain 

dv a I M-^m b 

—rr = V — Pb^ a W ab (12 

dt Pa ^ Pb 

for the acceleration of particle a. This form solves the Euler equation to the 
order of the method, but it docs not conserve the total momentum. Consider 
the force that particle b exerts on particle a 

( dv a \ m a m b 

Fba= [m a — ) = Pb^aW a b (13 

V dt J b Pa Pb 

and similarly, the force from particle a on b 

Fab= rnamb paVaW ^ 

Pa Pb 

where we have used Eq. (5). Since in general P a ^ Pf>, the sum over all the 
momentum derivatives, ^2 b d(mbV b )/dt, does not vanish and therefore the 
total momentum is not conserved. 

This deficiency can be easily cured by expressing the pressure gradient term 
as 
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^-4V,W?]. (15) 



p p \p 

If the gradient formula, Eq. (8), is applied to Eq. (15), the momentum equation 
reads 



t~s>(iHh"'* (16) 



Because the pressure part of the equations is now manifestly symmetric in a 
and b and V ' a W a b — — ^bW ba the forces are now equal and opposite ("actio= 
reactio" ) and therefore the total momentum is conserved by construction, i.e. 

Similarly, the total angular momentum is conserved since the sum of all 
torques vanishes: 



dt 

a,b 




" L = ^2r a x F ba = ^ \^r a x F ba + ^2r a x F 



ha 



l\52(r a -r b )xF ba )=0. (17) 



Here the summation indices were relabeled and F ab = —F ba was used. The 
expression finally vanishes, because the forces between particles act along the 
line joining them, see Eq. (4). 

A suitable energy equation can be constructed from Eq. (11) in a straight 
forward way. Start from the (adiabatic) first law of thermodynamics 



du a = Padfa 
dt p\ dt 



: u ( 18 ) 



and insert 

^2 m b W ab = ^ m b V ab ' V a VF a6 , (19) 




where we have used Eq. (6), to find 

du a _ _ 

n 2 • 

b 



n -^y^^bVab-VaWab. (20) 

dt p\ ^ 



Together with an equation of state the equations (7), (16) and (20) form a 
complete set of SPH equations. 

In the previous derivation it was implicitcly assumed that derivatives of the 
smoothing lengths can be ignored. In a simulation with strongly changing ge- 
ometry, however, it is advisable to locally adapt the smoothing length. This 
introduces, in principle, additional terms in the SPH equations. The impor- 
tance of these extra terms depends very much on the exact application [10, 17]. 
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The SPH-equations from a Lagrangian, "grad-h" terms 

The Lagrangian and the Euler-Lagrange equations 

The SPH equations can be derived by using nothing more than a suitable 
Lagrangian, the first law of thermodynamics and a prescription on how to 
obtain the density via summation. The Lagrangian of a perfect fluid [18] 



p^--u(p,s)jdV, (21) 

with s being the specific entropy, can be SPH-discretized in a straightforward 
way: 

£sPH,h = X] m6 (~2 - u (Pb, s b)j ■ (22) 



The discretized equations for the fluid are then found by applying the Euler- 
Lagrange equations 

d ( dL\ dL 



dt \dv a J dr a °' < ~ 23 '' 
The term in brackets yields the canonical particle momentum 

dL 

= m a v ai (24 

ov a 

the potential-type second term in the Lagrangian becomes 

dL v - du(p b , Sb) \ - dub 
dr a 2-j b Q ra 2-j b Q pb 



dpb 
dr a ' 



(25) 



The first derivative can be expressed using the first law of thermodynamics, 
du = P/ p 2 dp, and therefore 

dv a ^ P b dp b . . 



The density, its derivatives and the "grad-h" -terms 
We will now address the aditional terms resulting from variable smoothing 
lengths. For a density estimate as "local" as possible we use the smoothing 
length h a in 

Pa = ^2m b W{r ab ,h a ). (27) 

b 

Adaptivity can be reached by evolving the smoothing length according to 

h a = r,(^j'\ (28) 
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where r\ is a parameter typically in a range between 1.2 and 1.5[19] . Since p a 
and h a mutually depend on each other, see Eqs. (27) and (28) , an iteration is 
required for consistency. 

If we take the changes of h into account, the Lagrangian time derivative of 
the density is given by 

dp a ( dWab(ha) dr ab dW ab (h a ) dh a \ 



dt ^ mfc { dr ab dt dh a dt J 

= y^m b v ab ■WaW ab (h a ) + ^^^y^mb^-Wahtha), 
^ dp a dt ^ dh a 

where we have used dr ab /dt — e ab ■ v ab and Eq. (5). If the dp a /dt-terms are 
collected into the quantity 

^-(l-^-E-^^))' (29) 
the time derivative of the density reads 

^T = ^-J2m b v ab -V a W ab (h a ). (30) 

a b 

This is the generalization of the standard SPH expression, Eq. (19). 
In a similar way the spatial derivatives can be calculated 

dp b ^ f dW bk {h b )dh b \ 
_ = ^ mfe ^7 a W bk (h b ) + — ^ — | 

or, 

^ = ^E ro * V -^ fc (/i6). (31) 



+ 



The SPH equations with "grad-h" -terms 
Inserting Eq. (30) into Eq. (18) yields the "grad-h" energy equation 

= TT^Y,™^ ■ V a W ab (h a ). (32) 
dt S2 a p 2 a ^ 

With the derivative Eq. (31) one can write Eq. (26) as 

m a <h JT = - ^mh^VaPh = -E™ 6 ^! ( J)-^2 m k^aW bk (h b ) ] . (33) 

at b Pb b Pb ytb k j 
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With Eq. (3), the above equation becomes 
dv a Pb 1 



m a — f = -^2m b ^—^2m k V b W kb (h b ) (S ba - 5 ka ) 
at h p b ii b k 



= -m a ]T m b ( ^ V a W ab (h a ) + ^JaW ab (h b ) ) , (34) 



b 

i.e. the final momentum equation reads 
dv a>h 



dt 



E mb (-^^aWabiha) + J^\7 a W ab (h b )) . (35) 
h \ i2 aP a "bP b / 



Together with the density equation, Eq. (7), and an equation of state, Eqs. (32) 
and (35) form a complete set of "grad-h" SPH-equations. 

Self-gravity and gravitational softening 
The variational concept can also be applied to derive the gravitational forces 
including softening in a self-consistent way [20]. If gravity is taken into account, 
a gravitational part has to be added to the Lagrangian, Lsph = -ksPH.h + 
£sPH\ g with 

£sPH, g = -^2m b $b, (36) 

b 

where <l> b is the potential at the particle position b, <?(r b ). The potential # 
can be written as a sum over particle contributions 

^(r) = -G^m b <P(\r-r b \,h), (37) 

b 

and it is related to the matter density by Poisson's equation 

V 2 <£ = AirGp. (38) 

If we insert the sum representations of both the potential, Eq. (37) , and the 
density, Eq. (7), into the Poisson equation, Eq. (38), we obtain a relationship 
between the gravitational softening kernel, 4>, and the SPH-smoothing kernel 
W: 

W(\r - r b \,h) = (r 2 l r H\r - r 6 |, *)) • (39) 

Here we have used that both (f> and W depend only radially on the position 
coordinate. 

Applying the Euler-Lagrange equations, Eq. (23), to £ gr av,g yields the particle 
acceleration due to gravity [20] 
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dVa,g 

dt 



m b 



<P'ab( h a) + €b(hb) 



e a b 



G 



E 



m h 



^V a W ab (h a ) + ^-V a W ab (h b ) 
\l a s/ b 



(40) 



where <fi' ab — d<fi/d\r a — r b \. The first term in Eq. (40) is the gravitational 
force term usually used in SPH. The second term is due to gradients in the 
smoothing lengths and contains the quantities 



dhk d(j) k b{hk) 
Cfe = ^r}_^ m b- 



dp k 



dh k 



(41) 



and the flk defined in Eq. (29). Formally, it looks very similar to the pressure 
gradient terms in Eq. (35) with G£fc/2 corresponding to Pkj p\- As £fc is a 
negative definite quantity, these adaptive softening terms act against the gas 
pressure and therefore tend to increase the gravitational forces. The explicit 
forms of <j>, </>' and d(f)/dh for the cubic spline kernel an be found in Appendix 
A of [20]. 



2.2 Ideal magnetohydrodynamics 

Magnetic fields pervade the Universe in substantial strengths on all scales[21]. 
They are observed in intra-clustcr media in galaxy clusters [22] as well as in 
individual galaxies [23]. They are thought to be important for the birth of stars 
[24], they influence the life of stars e.g. via Sun spots or via controlling the 
angular momentum evolution during a stellar lifetime [25]. Stellar corpses such 
as neutron stars make themselves known via their magnetic field as pulsars, in 
a particular breed of neutron stars, so-called "magnetars" [26], the field reaches 
gigantic field strengths of the order <~ 10 15 Gauss. On the scale of planets, 
magnetic fields controle the magnctosphcrcs that can shield the planet from 
the lethal cosmic rays, a fact that has certainly facilitated the evolution of life 
on our planet. 



Basic equations of ideal MHD 

Magnetohydrodynamics is a one-fluid model for a highly conducting plasma. 
It assumes that electromagnetic fields are highly coupled to the electron-ion 
component so that if the fields have a typical frequency w and wave num- 
ber k, they fulfill uot^ ~ 1 and fcAh ~ 1, where Th and Ah are the typical 
hydrodynamic time and length scales. If 

A>ias (-0 < (m c ) (r h ) <Cl ' ^ 

where /3 p i as is the ratio between gas and magnetic pressure, r\ Ji the Larmor 
radius of the ions, m; and m e the ion and electron masses and n is the typical 
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ion collision time, is fulfilled, the plasma can be described by the equations of 
ideal magnetohydrodynamics[27]: 



dp 
dt 

dv l 1 dS ij 



P V ■ v (43) 

(44) 



dt p dxi 

du P— , 

= V • v 45 

dt p 

^ = -B(V-v) + (B-V)v, (46) 
where the magnetic stress tensor is given by 

S ij = _ PS ij + _L ^ B i B j _ ^B 2 S 1 ^ (47) 

and the B k are the components of the magnetic field strength. Note that for 
ideal magnetohydrodynamics only the momentum equation has to be modi- 
fied, the energy and the continuity equation are identical to the case of van- 
ishing magnetic field. The form of the momentum equation employed here 
formally accounts for _B(V • B) terms which are needed for momentum con- 
servation in shocks but on the other hand can be the cause of numerical 
instabilities, see [28] for a detailed discussion. 

Due to its relative simplicity in comparison to a more sophisticated plasma 
treatment magnetohydrodynamics and in particular ideal magnetohydrody- 
namics has been employed throughout a broad range of applications with 
sometimes not sufficient consideration about its range of applicability. Whether 
the conditions of applicabilty [29, 27] really hold needs to be checked for each 
problem individually. 



Euler potentials 

Being dissipationless the ideal MHD equations are conservative which leads 
to some important implications, the most powerful of which is probably the 
frozen flux theorem [30] which states that the magnetic field is carried around 
by the plasma. This kinematic effect is due to the evolution equation of the 
magnetic field, Eq. (46), and represents the conservation of magnetic flux 
through a fluid element. In reality, i.e. in the presence of dissipative terms, 
some slippage between the magnetic field and the plasma will occur. 
The idea that the magnetic field lines are carried around by the flow is closely 
related to the concepts of Euler potentials [31] which are sometimes also re- 
ferred to as Clebsch variables. For a review on Euler potentials we refer to 
[32, 33]. The basic idea is to present the magnetic field by two scalar variables, 
a and (3 such that 
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B = Vax V/3. (48) 
From this definition it is obvious that 

B ■ Va = = B ■ V/3, (49) 

in other words: a and (3 are constant along each field line and can therefore be 
used as field line labels. This is graphically represented in Fig. 2. The frozen 




Figure 2. The intersection of a plane of constant a with a plane of constant f3 
labels a magnetic field line. 



flux property of ideal MHD the simply translates into advecting a and (3 with 
a Lagrangian fluid clement: 

^=0 and ^=0. (50) 
dt dt K ' 



The Eulcr potentials naturally relate to the magnetic: vector potential via 

A = aVf3 + V£ (51) 

or 

A = -/3Va + VV>, (52) 

where £ and "0 are arbitrary smooth functions. It is straightforward to check 
that both of the above forms of the vector potential yield the magnetic field 
via 

V x A = VaxV(3 = B. (53) 

Thus, the \7 B = 0-constraint that is otherwise very hard to fulfill in a particle 
method [34, 19, 35] can be hard-wired into the numerical scheme by using the 
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advected quantities a and [3 to construct the magnetic field via Eq. (48). 
This approach has the additional case that, as long as the magnetic field 
is not strong enough to substantially influence the dynamics of the plasma, 
i.e. in the high-/3 p i as -case, the evolution of different initial field configurations 
can be explored by just re-processing an existing simulation with different 
initial values for a and (3. To find Euler potential pairs for a given magnetic 
field configuration is however usually a non-trivial task due to the non-linear 
nature of Eq. (48). The Euler potentials for a dipole field are known, but for 
more complicated field geometries its usually a challenge to find an analytical 
expression for the Euler potentials. There are however numerical procedures 
to find suitable pairs of Euler potentials, see e.g. [36, 37]. 
In two dimensions, a magnetic field can be represented by 

a = A z p = z, (54) 

where A z is the z-component of a vector potential. 

Limitations of the Euler potential approach 

While the Euler potential approach makes some otherwise rather challenging 
problems such as magnetic field advection (see below) a trivial task, they have 
their own difficulties and limitations. 

First, the Euler potentials for a given field configuration are not uniquely 
determined [33]. Assume, for example, that one particular set of Euler poten- 
tials, a\ and f}\, is known. Then for a second set that is a function of the 
known ones, a 2 = a 2 (ai,0i) and f3 2 = /^(cki, one finds 

„ fda 2 d(3 2 d(3 2 da 2 \ ^ a 

and therefore a 2 and (3 2 will also be a set of Euler potentials for the same 
field as long as the term in brackets is equal to unity. 

Second, by their very nature the Euler potentials are restricted to the purely 
non-disspative case and thus they are not immediately suited to treat the case 
of dissipative effects in a plasma. 

Third, there are restrictions with respect to the magnetic field geometries that 
can be represented by Euler potentials. It is, for example, impossible to repre- 
sent a linked ploidal and toroidal field. Nevertheless, as will be demonstrated 
in Sec. 3. 3, on a large set of standard MHD-test problems the Euler potential 
approach yields excellent results. 

From a numerical point of view they involve higher-order derivatives, sec 
Eqs. (44) and (48) which is usually numerically challenging. However, as will 
be shown below, this not necessarily has to degrade the accuracy of the solu- 
tion. 
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2.3 Dissipative terms 

In both hydrodynamics and magnetohydrodynamics we are interested in prin- 
ciple in the non-dissipative cases, see Sec. 2.1 and 2.2. The corresponding 
equations, however, allow for discontinuous shock solutions which need to be 
captured in order to allow for a physically correct and numerically stable 
solution. This can be done by cither making use of the analytical solution 
by locally solving a Ricmann-typc problem or by artificially spreading the 
discontinuities to a numerically resolvable width which means making them 
continuous. This latter artificial viscosity approach is most often used in the 
context of smooth particle hydrodynamics, although Riemann-solver-type ap- 
proaches also do exist [38, 39]. 

A careful design of artificial dissipation terms is essential to capture physically 
correct solutions. This was recently demonstrated at the example of Kclvin- 
Hclmholtz instabilities [40] . In the design of artificial dissipation terms we are 
guided by two principles: a) we want to use a form of the artificial dissipation 
equations that is oriented at Riemann-solvers [41] and b) we aim at applying 
dissipative terms only where they are necessary, i.e. near discontinuities, and 
follow in this respect Morris and Monaghan [42] who suggested to use time 
dependent dissipation parameters. 

Based on the analogy with Riemann solvers Monaghan [41] presented a gen- 
eral formulation of dissipative terms. It was noted that the evolution equations 
of every conservative quantity should contain dissipative terms to controlc 
discontinuities. This approach has been applied to ultra-relativistic [43] and 
magnetohydrodynamic shocks [35]. The "discontinuity capturing" term for a 
variable A is of the form 



where a a is a number of order unity that specifies the exact amount of dissipa- 
tion, v s i gt A is an appropriate signal velocity and p a b the average mass density 
of particles a and b. 

A comparison with the SPH expression for Laplacians [44] shows that the 
above equation is really an expression for [40] 



with rj oc a A v sig \r ab \. 

Following [42] the parameter that determines the exact values of the dissi- 
pative parameters, a a, is made time-dependent. This is put into effect by 
integrating an additional differential equation containing both a source term, 
Sa , that indicates the necessity of artificial dissipation and a decay term that 
contains the typical time scale, ta, it takes a particle to pass the discontinuity. 
The evolution equation of the dissipation parameter is given by 




(XAVsig.A 
Pab 




(56) 




(57) 
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daA,a CKA.a ~ OWn 

~T7~ = \-S A .a, (58) 

(it T A .a 

where a m i n is the minimum value to which we allow a a to decay. The decay 
time scale is given by 

TA,a = 7^— - — , (59) 

where C is a constant of order unity that is chosen after careful numerical 
experiments at problems with analytically known solutions. 



3 The MAGMA code 

Collisions between stars are very rare events in the solar neighbourhood. Close 
to centres of galaxies and globular clusters, however, the number densities of 
stars arc higher by up to a factor of 10 6 [45] and therefore stellar collisions are 
very common events. In fact, the innermost 0.3 lightyears of our Galaxy can 
be considered an efficient "stellar collider" [46] . A different type of encounter 
can occur for stellar binary systems that contain compact stellar objects. If 
born at close enough separations such systems can be driven towards merger 
by the emission of gravitional waves. Although rare per space volume these 
types of encounters release tremendous amounts of gravitational energy and 
are therefore potentially visibly out to cosmological distances thereby making 
huge volumes observationally accessible and producing a substantial observa- 
tional rate. 

Some of the most exciting astrophysical objects are thought to form in such 
encounters and since both neutron stars and white dwarfs are known to be 
threaded by very large magnetic fields, a careful study of such mergers requires 
the inclusion of magnetic fields and their evolution. 

3.1 Scope and physics modules 

The acronym MAGMA stands for a magnetohydrodynamics code for merger 
applications and this code has originally been developed for the study of mag- 
netized neutron stars [47, 48]. A very detailed description of this code can be 
found in [17]. 

For astrophysical studies the code contains several physics modules that go 
beyong the scope of this article and shall only be briefly sketched here. The 
interested reader is referred to the astrophysical literature. 

Equation of state 

For the thermodynamic properties of neutron star matter we use a temperature- 
dependent rclativistic mean-field equation of state [49, 50]. It can handle tem- 
peratures from to 100 MeV 6 , electron fractions from Y e ~ (pure neutron 



6 1 MeV corresponds to 1.16 • 10 10 K. 
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matter) up to 0.56 and densities from about 10 to more than 10 15 g cm~ 3 . 
No attempt is made to include matter constituents that are more exotic than 
neutrons and protons at high densities. For more details we refer to [51]. 

Neutrino emission 

The code contains a detailed multi-flavor neutrino leakage scheme. An addi- 
tional mesh is used to calculate the neutrino opacities that are needed for the 
neutrino emission rates at each particle position. The neutrino emission rates 
are used to account for the local cooling and the compositional changes due 
to weak interactions such as electron captures. A detailed description of the 
neutrino treatment can be found in [52]. 

Self-gravity 

The self-gravity of the fluid is treated in a Newtonian fashion. Both the grav- 
itational forces and the search for the particle neighbors are performed with 
a binary tree that is based on the one described in [53]. These tasks are the 
computationally most expensive part of the simulations and in practice they 
completely dominate the CPU-time usage. Forces emerging from the emission 
of gravitational waves are treated in a simple approximation. For more details, 
we refer to the literature [54, 51]. 



3.2 The MAGMA equations 

Here, we will only briefly summarize the implemented equations, the explicit 
forms of all the equations can be found in [17]. 

Instead of explicitely integrating the continuity equation, we calculate the 
density via summation as in Eq. (7). The momentum equation is used in the 
form 

dVgMRD _ cfaVh cfaq,h,diss dv q :g dv q^mag dv a, ma g,diss , fi „s 

dt _ dt dt dt dt dt [ ' 

where d(v a ^)/dt is given in Eq. (35), d(v a , g )/dt is given in Eq. (40), and the 
explicit forms of the dissipative terms, d(v a ^^i ss )/dt and c?(i>o,mag,diss)/cft, can 
be found in [17]. The magnetic force term is used in the form 



dt 



E — \§^^aW ab {h a ) + §4v a w ab (h b )\ 

//() ! I>nl>h } 



y, m b f B a (B a ■ V a W ab ) - B b (B b ■ V a W ab ) \ ^ 



(62) 



b 

where the symmetrized kernel gradient is given by 



V a W ab = 1 - 



^rV a W ab {h a ) + ^V a W ab {h b ) 
il a il b 
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The magnetic field is calculated from the Eulcr potentials 7 . Note that another 
form of the magnetic force term is also possible [35, 17]. The gradients of the 
Eulcr potentials are calculated in a way that gradients of linear functions are 
reproduced exactly [19, 17]. To handle magnetic shocks artificial dissipation 
terms were constructed according to the ideas outlined in Sec. 2.3. They are 
also applied to the evolution of a a and (3 a - They arc not meant to mimic 
physical dissipation in any way, their exclusive aim is to keep gradients nu- 
merically treatable. 

The MAGMA energy equation is of the form 

fifctq.MHD _ du aJl du a , AV du a>C 

dt dt + dt + dt ' ( ' 

where d(u a ,h) /dt is given in Eq. (32) , the explicit form of the artificial viscosity 
term d(u a Av)/dt and the thermal conductivity term d(u a c)/dt can be found 
in [17]. 

3.3 Tests and benchmarks 

We present here a selection of standard tests used in the hydro- and mag- 
netohydrodynamics community to validate numerical schemes. For a more 
exhaustive set of benchmarks we refer to [17] . 

Hydrodynamics 

ID: Sod's shock tube 

As a standard test of the shock capturing capability we show the results of 
Sod's shock tube test [55]. To the left of the origin, the initial state of the 
fluid is given by [p, P, uJl = [1.0,1.0,0.0] whilst to the right of the origin the 
initial state is [p,P,v x ]-n = [0.125,0.1,0.0] with 7 = 1.4. The problem is setup 
using 900 equal mass particles in one spatial dimension. Rather than adopting 
the usual practice of smoothing the initial conditions across the discontinuity, 
we follow [19] in using unsmoothed initial conditions but applying a small 
amount of artificial thermal conductivity. The results are shown in Figure 3, 
where the points represent the SPH particles. For comparison the exact solu- 
tion computed using a Riemann solver is given by the solid line. 
The shock itself is smoothed by the artificial viscosity term, which in this case 
can be seen to spread the discontinuity over about 6 particles. The contact 
discontinuity is smoothed by the application of artificial thermal conductiv- 
ity which (in particular) eliminates the "wall heating" effect often visible in 
numerical solutions to this problem. The exact distribution of particle sepa- 
rations in the contact discontinuity seen in Figure 3 is related to the initial 
particle placement across the discontinuity. 

7 Note that the code also allows to evolve magnetic fields according to a more 
straightforward SPH discretisation [35]. 
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Figure 3. Results of the Sod shock tube test in one dimension using 900 SPH 
particles setup using unsmoothed initial conditions. Artificial viscosity and thermal 
conductivity are applied to appropriately smooth the shock and contact discontinu- 
ity respectively. The exact solution is given by the solid line. The upper row displays 
the velocity (left) and the density (right), the bottom row shows specific internal 
energy (left) and the pressure (right). 



For this test, applying artificial viscosity and thermal conductivity as de- 
scribed, we do not find a large difference between the "grad-ft," formulation 
and other variants of SPH based on averages of the smoothing length. If any- 
thing, the "grad-ft"-terms tend to increase the order of the method, which, 
as in any higher order scheme, tends to enhance oscillations which may oth- 
erwise be damped, visible in Figure 3 as small "bumps" at the head of the 
rarefaction wave (in the absence of artificial viscosity these bumps appear as 
small but regular oscillations with a wavelength of a few particle spacings). 

3D: Sedov blast wave test 
In order to demonstrate that our scheme is capable of handling strong shocks 
in three dimensions, we have also tested the code on a Sedov blast wave prob- 
lem both with, see Sec. 3.3, and without magnetic fields. Without magnetic 
fields the explosion is spherically symmetric, however for a strong magnetic 
field the blast wave is significantly inhibited perpendicular to the magnetic 
field lines, resulting in a compression along one spatial dimension. Similar 
tests for both hydrodynamics and MHD have been used by many authors 
- for example by [56] in order to benchmark an Adaptive Mesh Refinement 
(AMR) code for MHD and by [57] in benchmarking the cosmological SPH 
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Figure 4. Results ol the hydrodynamic Sedov blast wave test in 3D at t = 0.09 
at resolutions of 125,000 (top) and 1 million (bottom) particles respectively. The 
density and radial position of each SPH particle are shown in each case, which may 
be compared to the exact solution given by the solid line. 



code GADGET. 

The hydrodynamic version is set up as follows: The particles are placed 
in a cubic lattice configuration in a three dimensional domain [—0.5,0.5] x 
[—0.5, 0.5] x [—0.5, 0.5] with uniform density p = 1 and zero pressure and tem- 
perature apart from a small region r < R near the origin, where we initialize 
the pressure using the total blast wave energy E = 1, ie. P = (7— 1)E/(^ttR 3 ). 
We set the initial blast radius to the size of a single particle's smoothing sphere 
R = 2rjAx (where 2 is the kernel radius, r)(= 1.5) is the smoothing length in 
units of the average particle spacing as in Eq. (28) and Ax is the initial parti- 
cle spacing on the cubic lattice) such that the explosion is as close to point-like 
as resolution allows. Boundaries are not important for this problem, however 
we use periodic boundary conditions to ensure that the particle distribution 
remains smooth at the edges of the domain. 

The results shown in Figure 4 at t = 0.09 have been obtained with a resolu- 
tion of 50 and 100 particles 3 (ie. 125,000 and 1 million particles respectively), 
where we have plotted (left panels) the density in a z — cross section slice 
and (right panels) the density and radial position of each particle (dots) to- 
gether with the exact self-similar Sedov solution (solid line). 
We found that the key to an accurate simulation of this problem in SPH is 
to incorporate an artificial thermal conductivity term due to the huge ini- 
tial discontinuity in thermal energy. The importance of such a term for shock 
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problems in SPH has been discussed recently by [19, 40]. In the absence of this 
term the particle distribution quickly becomes disordered around the shock 
front and the radial profile appears to be noisy. From Figure 4 we see that at 
a resolution of 1 million particles the highest density in the shock at t = 0.09 
is Pmax = 2.67 whereas for the lower resolution run p max = 2.1, consistent 
with a factor of 2 change in smoothing length. Using this we can estimate 
that a resolution of <~ 345 3 = 41 million particles is required to fully resolve 
the density jump in this problem in three dimensions. Note that the minimum 
density obtained in the post-shock rarefaction also decreases with resolution. 
Some small-amplitude post-shock oscillations are visible in the solution which 
we attribute to interaction of the spherical blast wave with particles in the 
surrounding medium initially placed on a regular (Cartesian) cubic lattice. 

Magnetohydrodynamics 

ID: Brio-Wu shock tube test 

The magnetic shock tube test of [58] has become a standard test case for nu- 
merical MHD schemes that has been widely used by many authors to bench- 
mark (mainly grid-based) MHD codes [59, 60, 61, 62]. The Brio-Wu shock 
test is the MHD analogon to Sod's shock tube problem that was described 
earlier, but here no analytical solution is known. The MHD Riemann problem 
allows for much more complex solutions than the hydrodynamic case which 
can occur because of the three different types of waves (i.e. slow, fast and 
Alfvcn, compared to just the sound waves in hydrodynamics). In the Brio-Wu 
shock test the solution contains the following components (from left to right 
in Fig. 5): a fast rarefaction fan and a slow compound wave consisting of a 
slow rarefaction attached to a slow shock (moving to the left) and a contact 
discontinuity, a slow shock and a fast rarefaction fan (moving to the right). 
It has been pointed out, however, that the stability of the unusual compound 
wave may be an artifact of the restriction of the symmetry to one spatial 
dimension whilst allowing the magnetic field to vary in two dimensions, [63] . 
The shown results are obtained using Eulcr potential formulation. Results 
of this problem using Smoothed Particle Magnetohydrodynamics (SPMHD) 
have been presented elsewhere [28, 19]. The Euler potentials show a dis- 
tinct improvement over the standard SPMHD results. The initial conditions 
on the left side of the discontinuity are [p,P,v x ,v y ,B y ] L = [1,1,0,0,1] and 
[p, P, v x , v y , B v ]r = [0.125,0.1,0,0,-1] on the right side. The x— component 
of the magnetic field is B x — 0.75 everywhere and a poly tropic exponent 
of 7 = 2.0 is used. Using the Euler potentials the components are given by 
a = —B y x (equivalent to the vector potential A z ) and (3 = z (or more specif- 
ically V/3 — z) and the B x component is treated as an external field which 
requires adding a source term to the evolution equation for a. Particles are 
restricted to move in one spatial dimension only, whilst the magnetic field is 
allowed to vary in two dimensions (that is, we compute a v y but do not use 
it to move the particles). 
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We setup the problem using 631 equal mass particles in the domain x G 
[—0.5,0.5] using, as in the hydrodynamic case, purely discontinuous initial 
conditions. Artificial viscosity, thermal conductivity and resistivity are ap- 
plied. The results are shown at t = 0.1 in Figure 5. For comparison the 
numerical solution from [62] is given by the solid line (no exact solution exists 
for this problem). The solution is generally well captured by our numerical 
scheme. Two small defects are worth noting. The first is that a small offset is 
visible in the thermal energy this is a result of the small non-conservation in- 
troduced by use of the Morris formulation [34] of the magnetic force, Eq. (61). 
Secondly, the rightmost discontinuity is somewhat over-smoothed by the arti- 
ficial resistivity term. We attribute this to the fact that the dissipative terms 
involve simply the maximum signal velocity v S i g (that is the maximum of all 
the wave types). Ideally each discontinuity should be smoothed taking ac- 
count of it's individual characteristic and corresponding v s i g (as would occur 
in a Godunov-MHD scheme). Increasing the total number of particles also 
decreases the smoothing applied to this wave. 




Figure 5. Results of the Brio & Wu MHD shock tube test at t = 0.1 using 631 
particles and the Euler potential formulation. For comparison the numerical solution 
taken from [62] is given by the solid line. The solution illustrates the complex shock 
structures which can be formed due to the different wave types in MHD, including 
in this case a compound wave consisting of a slow shock attached to a rarefaction 
wave. 



2D: Current loop advection problem 
A simple test problem for MHD is to compute the advection of a weak mag- 
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Figure 6. Magnetic field lines in the current loop advection test, plotted at t 
(top) and after 1000 crossings of the computational domain (bottom). 



netic field loop. This test, introduced by [64] in the development of the Athena 
MHD code 8 , presents a challenging problem for grid-based MHD schemes re- 
quiring careful formulation of the advection terms in the MHD equations. For 
our Lagrangian scheme, this test is straightforward to solve which strongly 
highlights the advantage of using a particle method for MHD in problems 
where there is significant motion with respect to a fixed reference frame. 
We setup the problem following [64] : the computational domain is two dimen- 
sional with x G [—1,1], y G [—0.5,0.5] using periodic boundary conditions. 
Density and pressure arc uniform with p = 1 and P = 1. The particles are 
laid down in a cubic lattice configuration with velocity initialized according 
to v = (v cos 9, v sin 9) with cos# = 2/v / 5, sin# = l/\/5 and v — 1 such 
that by t = 1 the field loop will have been advected around the computational 
domain once. The magnetic field is two dimensional, initialized using a vector 
potential given by 

= (A (R-r)r<R, 
Az -\0 r>R, (M) 

where Aq = 10~ 3 , R = 0.3 and r = y 1 x 1 + y 2 . The ratio of thermal to mag- 
netic pressure is thus given by f3 p \ as = P/(\B 2 ) = 2 x 10 6 (for r < R) such 
that the magnetic field is passively advected. [64] show the results of this 
problem after two crossings of the computational domain, by which time the 
loop has either been significantly diffused or has disintegrated into oscillations 

8 http://www.astro.princeton.edu/~jstone/athena.html 
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depending on details of their particular choice of scheme. The advantages of 
a Lagrangian scheme are that advection is computed exactly, and using our 
Eulcr potential formulation for the magnetic field (which in two dimensions 
is equivalent to a vector potential formulation with a — A z and (3 = z), this 
is also true for the evolution of the magnetic field. The result is that the field 
loop is advected without change by our code for as long as one may care to 
compute it. This is demonstrated in Fig. 6 which shows the magnetic field 
lines at t = (top) and after 1000 (!) crossings of the computational domain 
(bottom) , in which the field configuration can be seen to be identical to the top 
figure. The magnetic energy (not shown) is also maintained exactly, whereas 
[64] find of order a 10% reduction in magnetic energy after two crossings of 
the domain. 

In a realistic simulation involving MHD shocks there will be some diffusion 
of the magnetic field introduced by the addition of artificial diffusion terms, 
which are required to resolve discontinuities in the magnetic field. However 
the point is that these terms are explicitly added to the SPH calculation and 
can be turned off where they are not necessary whereas the diffusion present 
in a grid-based code is intrinsic and always present. 

2D: Orszag-Tang test 
The evolution of the compressible Orszag-Tang vortex system [65] involves 
the interaction of several shock waves traveling at different speeds. Originally 
studied in the context of incompressible MHD turbulence, it has later been 
extended to the compressible case [66, 67]. It is generally considered a good 
test to validate the robustness of numerical MHD schemes. In the SPH con- 
text, this test has been discussed in detail by [19] and [35]. 
The problem is two dimensional with periodic boundary conditions on the do- 
main [0, 1] x [0, 1]. The setup consists of an initially uniform state perturbed 
by periodic vortices in the velocity field, which, combined with a doubly peri- 
odic field geometry, results in a complex interaction between the shocks and 
the magnetic field. 

The velocity field is given by v = t>o[— sin (27ry), sin (27ra)] where vo = 1. The 
magnetic field is given by B — B [— sin (2ny), sin (4ttx)} where B = 1/VAn. 
Using the Euler potentials this corresponds to a = A z = £> /(27r)[cos (27ry) + 
| cos (47rx)]. The flow has an initial average Mach number of unity, a ra- 
tio of magnetic to thermal pressure of 10/3 and we use a polytropic expo- 
nent 7 = 5/3. The initial gas state is therefore P = 5/3Bq = 5/(127r) and 
p = jP/vq = 25/(367r). Note that the choice of length and time scales differs 
slightly between various implementations in the literature. The setup used 
above follows that of [61] and [68]. 

We compute the problem using 512 x 590 particles initially placed on a uni- 
form, close-packed lattice. The density at t = 0.5 is shown in Figure 7 using 
both the SPMHD formalism of [35] (left), and the Euler potential approach 
(right) outlined in Sec. 2. 2. The Euler potential formulation is clearly supe- 
rior to the standard SPMHD method. This is largely a result of the relative 
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Figure 7. Density distribution in the two dimensional Orzsag-Tang vortex prob- 
lem at t = 0.5. The initial vortices in the velocity field combined with a doubly 
periodic field geometry lead to a complex interaction between propagating shocks 
and the magnetic field. Results are shown using 512 x 590 particles using a SPMHD 
formalism of [35] (left) and using the Euler potentials (right). The reduced artifi- 
cial resistivity required in the Euler potential formalism leads to a much improved 
effective resolution. 



requirements for artificial resistivity in each case. In the standard SPMHD 
method the application of artificial resistivity is crucial for this problem (that 
is, in the absence of artificial resistivity the density and magnetic field distri- 
butions are significantly in error). Using the Euler potentials we find that the 
solution can be computed using zero artificial resistivity, relying only on the 
"implicit smoothing" present in the computation of the magnetic field using 
SPH operators. This means that topological features in the magnetic field 
are much better preserved, which is reflected in the density distribution. For 
example the filament near the center of the figure is well resolved using the 
Euler potentials but completely washed out by the artificial resistivity in the 
standard SPMHD formalism. Also the high density features near the top and 
bottom of the figure (coincident to a reversal in the magnetic field) arc much 
better resolved using the Euler potentials. 

3D: MHD blast wave 
The MHD version of the Sedov test is identical to the hydrodynamic test 
with the addition of a uniform magnetic field in the a;— direction, that is B = 
[Bo, 0,0] with Bq = 3.0. Initially the surrounding material has zero thermal 
pressure, meaning that the plasma /3 p i as is zero (ie. magnetic pressure infinitely 
strong compared to thermal pressure). However, this choice of field strength 
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Figure 8. Results of the 3D MHD blast wave test at t = 0.05 at a resolution of 1 
million (100 3 ) particles. Plots show (left to right, top to bottom) density, pressure, 
magnitude of velocity and magnetic field strength (with overlaid field lines) , plotted 
in a cross-section slice through the z — plane. 



gives a mean plasma /3 p i as in the post-shock material of /3 p i as ~ 1.3, such that 
the magnetic pressure plays an equal or dominant role in the evolution of the 
shock. The results of this problem at t = 0.05 are shown in Fig. 8, where plots 
show density, pressure, magnitude of velocity and magnetic field strength in 
a cross section slice taken at z = 0. In addition the magnetic field lines are 
plotted on the magnetic field strength plot. 

In this strong-field regime, the magnetic field lines arc not significantly bent by 
the propagating blast wave but rather strongly constrain the blast wave into 
an oblate spheroidal shape. The density (and likewise pressure) enhancement 
in the shock is significantly reduced in the y— direction (left and top right 
panels) due to the additional pressure provided by the magnetic field which 
is compressed in this direction (bottom right panel). 



4 Summary and conclusion 

We have outlined several recent developments in smooth particle hydrody- 
namics. The equations of self-gravitating, ideal hydrodynamics were derived 
cxplicitely from a Lagrangian thereby yielding the correct particle index sym- 
metries that ensure that the physical conservation laws are hard-wired into 
the discrete set of SPH equations without any arbitrariness. We have further 
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described the implementation of ideal MHD via so-called Euler potentials. 
This approach enforces the crucial V • B = O-constraint by construction. 
All dissipative terms required to capture discontinuities were carefully de- 
signed so that they a) have a form suggested in analogy with Riemann-solvers 
and b) are only active near discontinuities. These principles are implemented 
in our three-dimensional, Lagrangian magnetohydrodynamics code MAGMA. 
In a set of standard test problems used to benchmark numerical (magneto- 
hydrodynamics schemes we have demonstrated the excellent performance of 
the code. 
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